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INTERACTIVE  FLUTTER  SOLVER  AT  ARL 


by 

S.  GALEA 


SUMARY 

An  interactive  direct  flutter  solving  routine  has  been  Installed  on 
ARL's  ELXSI  computer.  Interactive  graphics  routines,  using  DI-3000 
graphics,  software,  have  been  incorporated  to  give  the  user  a 
progressive  picture  of  the  solution.  Estimates  of  subcrltlcal  response 
data  may  also  be  determined  using  this  method.  A  description  of  the 
method  and  the  associated  software  is  presented  here.  Also  included  is 
a  two-dimensional  flutter  problem,  using  quasi-steady  aerodynamics, 
which  is  solved  by  this  direct  flutter  solver. 
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1. 


1 .  INTRODUCTION 

The  classical  methods  for  the  solution  of  the  flutter  equations  are 
the  k  method  ("American  Method")  and  the  p-k  method  ("British  Method"), 
see  [1,2].  Both  these  methods  treat  the  flutter  equations,  which  are  a 
set  of  n  second-order  linear  homogeneous  differential  equations,  as  an 
eigenvalue  problem  and  solve  for  the  eigenvalue,  x*  Here  the  imaginary 
component  of  x  gives  the  frequency  of  the  response  while  the  real 
component  provides  the  decay  rate.  Flutter  occurs  when  the  oscillation  is 
non-convergent ,  i.e.,  when  the  decay  rate  is  zero. 

The  flight  speed  at  which  this  occurs  is  called  the  flutter  speed. 
The  k  method  assumes  that  the  aircrafts  aeroelastic  response  is 

sinusoidal  and  hence  the  solution  is  correct  only  at  the  flight  flutter 
speed.  No  such  restriction  is  made  in  the  p-k  method.  The  latter  method 
does  assume,  however,  that  the  aerodynamic  forces  due  to  constant 

amplitude  oscillating  lifting  surfaces  are  equal  to  aerodynamic  forces 
generated  by  lifting  surfaces  with  slowly  increasing  or  decreasing 
amplitudes  of  oscillation.  This  assumption  is  accurate  for  low  values  of 
decay  rate,  either  positive  or  negative.  The  p-k  method  is  therefore  more 
accurate  in  the  subcrltlcal  and  super  critical  regime  than  the  k  method. 

In  flutter  problems  the  prime  concern  is  to  calculate  the  flutter 

speed,  in  which  case  the  k  and  the  p-k  methods  will  provide  accurate 
results.  Another  concern  is  the  aircraft's  suboritical  flight 
characteristics,  whereby  if  in  flight  or  wind  tunnel  tests  the  flutter 

speed  is  not  attained  then  these  characteristics  can  be  used  for 
comparison  with  the  appropriate  mathematical  model.  In  this  case  the  p-k 
method  will  give  more  accurate  results  than  the  k  method. 

A  new  method  of  solving  the  flutter  equation  is  described  in 
References  3  and  4.  The  present  paper  describes  a  computer  program  based 
on  this  method  which  provides  estimates  of  subcrltlcal  response  data. 
Unlike  classical  methods,  this  direct  flutter  solver  allows  the  user  to 
establish  quickly  relationships  between  any  two  parameters  which  satisfy 
the  flutter  equations. 


2. 


2.  THE  FLUTTER  EQUATION 

The  aeroelastic  motion  of  an  aircraft  can  be  modelled  by  the  flatter 
equations  which  may  be  given  in  non-dimensional  form  as: 

[A]  q+[B  +  D]q+[C  +  E]q»0  (1) 

where  A,  D,  and  E  are,  respectively,  the  structural  inertia,  damping  and 
stiffness  matrices,  B  and  C  are,  respectively,  the  aerodynamic  damping  and 
stiffness  matrices,  and  q  is  the  generalised  co-ordinate  vector.  The 
aerodynamic  matrices  are  functions  of  the  non-dimensional ised  frequency  v, 
Mach  number  M,  velocity  U,  air  density  p  ,  and  aircraft  planform. 

t,  b)& 

Here  v  ”  y“ 

where  u  is  the  frequency  of  oscillation,  1  is  the  reference  length  and  U 

the  flight  velocity.  In  this  case  q  -  Bq/dt  where  -t  is  the  non- 

dimensional  time  given  by  -r  -  tU/t.  On  substituting  q  (t)  •  qQ  and 
dividing  by  equation  1  becomes  - 

1  X*  CA3  ♦  X  C  B  ♦  D  ]  ♦  [  C  ♦  E  ]  }  qQ  -  0  (2) 

Here  x  ■  6  *  iv,  6  -  vA  and  A  is  the  decay  rate  of  the  response. 

Equation  2  may  be  rewritten  as  - 

[S]  q0  -  0  (3) 

where  S,  the  system  matrix,  is  a  function  of  U,  M,  v,  air  density  (p), 
stiffness,  damping,  mass  distribution  and  so  on. 

3.  METHOD  OF  SOLUTION 

3.1  Refining  Initial  Estimatea 

Consider  n  real  parameters  x  -  x2,  ....  xR  which  specify  the 

model t  then  eqn  (3)  oan  be  expressed' as i 


(*) 


[  S(x)  ]  qe  ■  0 
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For  a  non-trlvlal  solution  the  aystea  aatrlx  must  be  singular;  thus  the 
determinant,  which  is  normally  complex,  must  be  zero. 

Following  the  notation  in  Reference  [4],  if  two  of  the  system's 
parameters,  Xj  and  xk,  are  not  fixed  then  the  equation  that  must  be 
satisfied  is 

|  Six)  |  -  F(Xj,xk)  -  F  <Xj,xh)  ♦  i  F  <XjiXh)  -  0  (5) 

Gaussian  elimination  is  used  to  evaluate  the  determinant  of  the 
system  matrix  S.  For  large  degree  of  freedom  systems  the  slowest  step  in 
evaluating  FtxjX^)  is  the  evaluation  of  the  system  matrix  S. 

To  minimize  the  number  of  function  evaluations,  the  secant  method  is 
used  to  solve  for  the  two  variables, Xj  and  x^,  so  as  to  satisfy  equation 
(5).  This  is  described  more  fully  in  Reference  4. 

The  secant  method  requires  three  points  at  any  time  to  evaluate  a 
refined  estimate.  For  the  first  iteration,  the  three  points  are  given 
by  F(xj  ,x*  ),  F(Xj  ♦  &Xy  x*)  and  F(xj,  x*  ♦  x*  ♦  6xk)  in  the 

F,  Xj,  xk  co-ordinate  system  and  F  (xj,  x£),  F*  (xj  ♦  6Xj,  x*)  and 
“F  (xj,  x*  ♦  «xk)  in  the  "f,  Xj,  xk  co-ordldnate  system.  The  initiate 
estimate  is  denoted  by  the  superscript  zero.  These  points  then  determine 
a  plane  in  each  co-ordinate  system.  The  refined  estimate  (xj,  x£)  is 
obtained  by  determining  the  point  of  intersection  of  the  above  two  planes 
with  the  plane  F  -  f  •  0. 

The  solution  Is  said  to  have  converged  at  the  nth  iteration  when  - 


lo«10C( • 


I  ♦  I 


k  ,--k  )  X  100  3  <  0.0002 

xn 

*k 


(6) 


If  convergence  does  not  occur  at  the  nth  iteration  then  one  of  the  points 

(x"~3.  x"”3),  (*r\  x”“2)  or  (X*?'1.  x”“ 1 ) ,  whichever  is  furthest  from 
J  K  v  K  n  n  J  K 

the  refined  estimate  (x. ,  x.  ),  is  discarded  and  so  the  two  surviving 

ft  ft  ^ 

points  plus  (x. ,  x‘)  are  used  to  obtain  a  further  refined  estimate 

na  1  J  * 
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3.2  Scaling 

In  evaluating  the  determinant  of  S,  computational  problems  may  arise 
either  when  its  diagonal  elements  are  very  large  leading  to  overflow,  or 
when  they  are  very  small  so  that  significant  figures  may  be  lost*  In 
order  to  overcome  these  problems  each  row  of  the  determinant  is  divided  by 
its  diagonal  element. 

3.3  Predicting  Initial  Estimates 

The  idea  of  this  direct  method  of  solution  is  to  solve  for  two 

dependent  parameters  Xj  and  xk  as  one  varies  a  particular  independent 

parameter  Xj .  That  is,  to  observe  the  effect  on  Xj  and  xk  as  x^  varies. 

Assume  a  particular  independent  parameter,  denoted  by  x<  has  the 

*  *p 

solution  (xj  p,  xk^p).  The  procedure  for  predicting  the  new  initial 

estimate  is  shown  in  Figure  1.  To  march  along  the  xi  axis,  a  small 

increment  is  closer  such  that  x,  .  »  x.  ♦  6x, .  Thus  the  initial 

estimate  (x^  xk  is  taken  as  (Xj>p,  xk,p^*  causing  rapid 

convergence  to  the  solution  xk,p«*1^#  Tlie  estimate  for  new 

point  xifp+2  (xj  x£  p+2^  13  obtained  linear  extrapolation 

using  the  previous  two  solutions. 

The  step  6x^  is  small  to  ensure  rapid  convergence  Once  convergence 

has  been  attained  an  arc  of  a  circle  is  fitted  through  the  three  points 

xi,p»  xi,p+1*  and  xi,p+2f  t0  ot>tain  an  estimate  for  the  dependent 
parameters  at  x<  v  Fitting  an  arc  of  a  circle  allows  one  to  rise 
large  6x^  and  still  obtain  aceurate  initial  estimates  (Xj,  x^)  which  are 
essential  to  obtain  rapid  convergence  to  the  solution  (x^x^). 

4.  PROGRAM 

Appendix  A  contains  descriptions  of  the  nain  subroutines  used  to 
solve  the  flutter  equation.  The  user-supplied  main  program  PREPFLSO  and 
subroutine  SYSCAL  are  described  in  Appendix  B.  A  block  diagram  of  all 
subroutines  used  in  program  VAFLUTSOLV  ia  shown  in  Figure  2. 


5.  EXAMPLE  TV.,_ 

m  considered  here.  Tnis 
a  al.ple  two-degree-of-fr.odc.  sya  validate  the 
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(Ta) 

.  n  *  sa ;  •  gh » * u  • 0 

am 

and 

_  .  mass  of  section, 

Her®  © 

K  .  bending  stiffness, 

^  «  torsional  stiffness, 

.  resultant  lift  force  per  unit  span 
r  .  moment  of  inertia  of  the  section  about  the  e.a.. 

*  -  moment  about  elastic  axis  due  to  L. 

J  .  distance  from  e.a.  to  aerodynamic  centre, 

g  .  mx  »  and 

x°  -  distance  from  e.a  to  c.g. 

°  ,  the  resultant  section  lift 

Using  quasi-steady  aerodynamic  the  y 

force,  L,  is  given  by 


q  c  l  a  ♦  -  3 

9a  « 


where 


1/2  p  U 

density  of  air, 
flight  velocity# 
wing  chord,  and 
lift  coefficient* 


6, 


*  .a 


Using  I  -  ac  r  ,  S  -  ■  x  c  and  Mv  -  Le  and  substituting  eqn  (8) 

Q  01  Cl  ft  J 

into  eqn  (7),  eqn  (7)  nay  be  written  as 


and 


-  t  »  h 

h  ♦  Ma  c  b  ♦  Kh  h  ♦  j  pU  c  — -[  «♦  j  3  «  0 

da 


a  _  1 


ac, 


c  r  a  ♦  m  x  c  h  ♦  K  a  -  5  pu  ce  -  C  a  ♦  s  ]  -  0 

a  a  0  2  3a  U 


(9) 


Let  U 
h 

T 
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a 

“h 
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*  .ft 

(U  /c)  h 

a  _n 
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K  /I 
a  a 


a  a  _> 

and  by  aultlplying  eqn  (9)  by  (1/2  pU  c  )  ,  eqn  (9)  can  be  non- 

dimen^ionalised  to  give,  in  aatrix  fora, 
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Note  that  the  overbar  denotes  non-diaensional  quantities. 
Asauaing  the  solution  is  oscillatory,  then 


and 


h  -  h0  l 
®  -  ®o  1 


Xt 


Xt 


(ID 
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where  x 
v 
A 


(Av  ♦  iv) 
wc/U 

decay  rate  (fraction  of  critical) 


Alternatively, 


X 


1 

0 


ie  x  X^U 


Using  eqns  (11)  and 


(12)  and  rearranging,  eqn  (10)  becomes, 


(12) 


where 


2m 

wpc 


The  above  eqn  is  of  the  form 


C  S(x)  3  qQ  -  0 


<1*0 


and  so  can  be  solved  by  our  method. 

Appendix  C  gives  a  list  of  the  parameters,  their  corresponding  values 

and  listings  of  the  user-supplied  routines.  Figure  4  contains  plots  of 

m/m  and  Am/m  versus  U  •  Comparisons  with  results  published  in  Dowell 
a  a 

et  al.  (Reference  [5])  are  also  shown  In  Figure  M,  and  It  can  be  seen  that 
a  good  agreement  has  been  obtained. 


6.  DISCUSSION 

The  following  Is  a  list  of  improremants  that  oould  be  incorporated 
into  the  preaent  program. 


8. 


1.  A  reasonable  initial  estimate  is  required  so  that  the  solution 
will  converge  quickly.  Consequently  one  could  incorporate  a 
classical  eigenvalue  solver  to  obtain  initial  solutions  of 
damping  and  frequency  at  V  -  0. 

2.  An  interpolation  procedure  to  determine  the  flutter  speed 
accurately. 

3.  Use  of  a  conic  or  cubic  spline  extrapolation  to  estimate  the  next 
Xj  and  xk.  Presently  only  an  arc  of  a  circle  is  used. 

There  are  a  number  of  options  contained  in  the  original  program 
(described  in  Reference  [H])  that  have,  as  of  yet,  not  been  included  in 
the  version  at  ARL.  These  options  are:- 

1.  Display  of  error  function  (see  eqn  (6))  against  number  of 
iterations. 

2.  Stability  test,  i.e.,  if  for  a  particular  solution  of  Xj  Vs  xif 
say,  the  domain  of  stability  is  not  particularly  clear.  This 
test  determines  which  domain  is  stable. 

7.  CONCLUSION 

A  direct  method  of  solving  the  flutter  equations  has  been  programmed 
at  ARL.  The  program  is  interactive,  Graphics  routines,  using  DI-3000 
software,  have  been  incorporated  to  give  the  user  an  immediate  picture  of 
the  solution's  progress.  The  user  may  evaluate  the  relationship  between 
any  two  parameters  for  which  the  flutter  equation  is  solved,  the  only 
prerequisite  being  that  these  are  variables  in  the  subroutine  that 
calculates  the  system  matrix  (i.e.  subroutine  SYSCAL).  This  method  also 
provides  estimates  of  subcritical  response  data. 
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APPENDIX  A 


Flutter  Solver  Routines 


Table  1  gives  a  list  of  the  routines  used  in  VAFLUTSOLV.  A  more 
detailed  discussion  of  each  routine,  except  for  the  user  supplied  routines 
(see  Appendix  B),  is  given  below. 

Subroutine  IN IT 

This  subroutine  prompts  the  user  on  - 

(a)  the  choice  of  independent  and  dependent  parameters  to  be  used 
initially,  and 

(b)  the  initial  estimates  of  the  dependent  parameters  for  a 
particular  independent  parameter. 

If  convergence  has  not  been  achieved  in  20  iterations  (this  check  is 
carried  out  in  subroutine  F2V)  then  the  user  has  the  option  of  entering 
another  estimate  or  continuing  with  the  present  estimate. 

Subroutine  OUTP 

This  subroutine  outputs  the  present  values  of  the  independent  and 
dependent  parameters. 

Subroutine  CHAMP  (OPTIONS,  IOPT,  ICON) 

When  called,  this  subroutine  displays  a  list  of  options  contained  in 
the  character  array  OPTIONS.  Generally  the  integer  variable  IOPT 
specifies  the  number  of  options  available  and  ICOM  is  the  integer  returned 
to  the  caller  routine  (i.e.  ICOM  is  the  option  selected).  Other 
variations  are,  if  on  entry 


IOPT 

IOPT 


0  then  the  DI  3000  graphics  software  is  initialized 
-1  then  the  screen  is  cleared,  and 


(A. 2) 


IOPT  -  -2  then  the  screen  is  set  up  for  interactive  plotting. 

If  on  calling  CMAND,  ICOM  -  -2,  then  options  are  set  up  but  no  command  is 
returned  to  the  caller  routine 

Subroutine  TPTBXT 

This  routine  provides  a  list  of  parameters  with  corresponding  values. 

Subroutine  PLABEL  (XMIN,  YMIN,  XMAX,  YMAX,  LABX ,  LABY,  SIZE) 

When  called,  this  subroutine  selects  suitable  scaling  for  the  x  and  y 
axes.  The  graph  axes,  with  appropriate  labelling,  are  drawn  and  stored  in 
the  integer  array,  PICT1 *  using  routines  from  GRAFMAKER  graphics  software 
(which  is  based  on  DI  3000  graphics  software). 

The  first  four  dummy  arguments,  listed  above,  give  the  lower  and 
upper  limits  of  the  variables  to  be  plotted  (here  x  and  y  denote  the 
horizontal  and  vertical  axes,  respectively)  Character  variables  LABX  and 
LABY  contain  the  x  and  y  label  names,  respectively.  The  desired  maximum 
dimension  of  the  plot  is  given  by  the  real  variable  SIZE,  where  the  screen 
is  assumed  to  be  10  units  square  (i.e.  SIZE  <  10). 

Subroutine  DISP 

On  calling  subroutine  DISP  a  plot  of  the  first  dependent  parameter 
versus  the  independent  parameter  is  displayed.  Subroutine  CMAND  is  then 
called;  thus,  the  following  menu  is  displayed: 

VRBLS  1 

SELECT  2 

C0MT  3 

DELETE  4 

STABLITY  5 

STEP  SZ  6 

SKP  PT  7 

PLOT  8 

CHG  VBL  9 


(A. 3) 


Once  an  option  is  selected,  subroutine  DISP  then  executes  the 
command.  A  description  of  each  option  is  given  below:- 

(1)  VRBLS  -  step  (a)  Displays  a  plot  of  the  second  dependent  parameter 

versus  the  independent  parameter, 

step  (b)  Displays  the  menu 

SHW  IVZ  1 
CHG  VBL  2 
CONT  3 

where  SHW  IVZ  again  displays  a  plot  of  the  first  dependent  parameter 
versus  the  independent  parameter,  CHG  VBL  allows  the  user  to  select  new 
Independent  and  dependent  parameters  and  CONT  causes  the  main  menu  to  be 
displayed  again. 


(2)  SELECT  Allows  the  user  to  select  new  starting  values.  This  is 

to  be  used  if  the  solution  has  not  converged  or  if  the 
dependent  and  independent  parameters  have  been  changed. 


(3)  CONT 


Returns  the  user  to  the  caller  routine  (DRIVNF), 


(4)  DELETE 


Deletes  the  previous  entry. 


(5)  STABLITY  Not  operative. 


(6)  STEP  SZ  Displays  the  menu, 


INCRSE  1 
DECRSE  -1 
SAME  0 


By  entering  1,  -1,  or  0  the  user  may  double,  halve  or  not  alter  the 
step  size,  respectively. 


(7)  SKPPT 


Not  operative 


(A  A) 


(8)  PLOT  Stores  plots,  presently  displayed  on  the  screen,  into  a 

metafile  (Refer  to  the  DI-3000  manual).  The  name  of  the 
metafile  is  specified  by  the  user  after  the  screen 
displays  the  prompt: 

...ENTER  FILE  NAME. 

(9)  CHG  VBL  Allows  the  user  to  select  new  independent  and  dependent 

parameters . 

Subroutine  TORT  (XO,  YO,  XI,  Y1 ,  X2,  Y2,  ICL) 

Given  co-ordinates  of  a  circle  centre  (XO,  YO)  and  two  points  (XI , 
Y1 ) ,  (X2,  Y2)  on  the  circle,  subroutine  VORT  determines  the  direction  from 
1  to  2.  Subroutine  VORT  returns  ICL  -  1  if  the  direction  is  clockwise  or 
ICL  -  -1  if  it  is  anticlockwise. 

Subroutine  AMO  (X,  Y,  TH) 

Given  the  co-ordinates  (X,  Y),  subroutine  ANG  calculates  the  angle  in 
radians  where  0  <  TH  <  2ir  . 

Subroutine  CIRCLE  (X,  Y,  XC,  YC,  RAD) 

This  subroutine  establishes  the  centre  of  a  circle  (XC,  YC)  that 
passes  through  three  points  contained  in  the  dummy  arrays  X(3)  and  Y(3). 
The  radius  of  the  circle,  RAD,  is  also  found. 

Subroutine  Iff 

This  subroutine  factorises  the  dependent  parameters,  using  the  real 
array  FACTR  (see  PREPFLSO),  and  then  calls  subroutine  F2V. 

Subroutine  FG  (XX,  Y,  F,  G) 

Given  dependent  parameters  XX  and  Y  this  subroutine  calls  subroutine 
SYSCAL  to  calculate  the  system  matrix  S.  Subroutine  CMINV  is  then  called 
to  obtain  the  complex  determinant  of  S  (l.e.  |s|  -  F  ♦  i  G). 
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Subroutine  CHIMV  (At  N,  ND,  B,  L,  0 ,  IRROR) 

If  L  is  negative,  zero  or  positive  then  the  N  x  N  matrix  A  is 
Inverted,  by  Gauss-Jordon  Pivotal  elimination,  and  A  is  replaced  by  its 
inverse. 

If  L  is  positive  then  the  N  x  L  matrix  B  is  manipulated  by  the  same 
procedure  as  used  to  invert  A  and  is  then  replaced  by  the  resulting 
matrix. 

Subroutine  F2V  (XX,  YY,  EPS,  NC,  IR,  FG) 

This  subroutine  calls  subroutine  FG  to  determine  the  determinant  |S| 
given  the  two  dependent  parameter  estimates  (XX,  YY).  These  estimates  are 
then  refined  by  using  the  secant  method.  Convergence  of  the  estimates  is 
checked  by  using  eqn  (6).  If  either  convergence  is  attained  or  the 
maximum  number  of  iterations  is  exceeded  then  the  user  is  returned  to  NF 
and  hence  to  DRIVNF. 

Subroutine  DRIVNF 

Subroutine  DRIVNF  is  the  main  subroutine  of  program  VAFLUTSOLV.  This 

subroutine  calls  INIT  to  obtain  the  initial  estimates  from  the  user.  That 

is,  for  a  particular  independent  parameter  value  (x<  0),  estimates  x  ° 

1 j  ,p 

and  x.°  are  supplied  by  the  user.  Subroutine  NF  is  then  called  (at 

K  ,  p 

least  three  times)  to  evaluate  solutions  for  the  first  two  Increments  in 
the  independent  parameter  (i.e.  solutions  are  obtained  at  x^p,  xi  p^i» 
and  Xj  in  a  manner  similar  to  that  discussed  in  Section  3* 

Subroutine  CIRCLE  and  VORT  are  then  called  to  determine  the  initial 
estimates  for  the  next  point  (xi#p+4)*  Once  again  NF  is  called  to  solve 
for  the  dependent  parameters  (xj#p«.4*  xk,p«V#  the  solution  haa  failed 
to  converge  (see  Subroutine  F2V)  then  the  step  size  of  the  independent 
parameter  is  halved,  and  the  above  procedure  is  repeated.  If  the  initial 
estimate,  at  half  the  step  size,  does  not  converge  then  the  step  size  is 
halved  again.  Rapidly  varying  functions  of  x^  usually  cause  estimates, 
when  fitting  an  arc  of  a  circle,  to  be  Inaccurate  and  hence  convergence  is 
hard  to  attain. 

If  convergence  has  not  occurred  after  carrying  out  the  above 
procedure  then  the  program  steps  backward,  in  small  steps,  to  obtain  two 
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other  solutions  just  upstream  of  the  last  convergent  point  (i.e.  in 

this  case).  With  these  three  points  a  new  estimate  (xj®p+ij*  xi<*ptV  ls 
evaluated,  again  fitting  an  arc  of  a  circle,  and  the  program  proceeds  as 
usual . 

Interactive  graphics  displays  are  obtained  by  calling  Subroutine 

DISP. 


APPENDIX  B 


User  Supplied  Routines 


Program  PBEPPLSO 

COMMON/VR/  VVRBL,  FACTR,  IVRBL,  IV1 ,  IV2,  NM 

COMMON/LAB  VRBL 


CHARACTER  *10  VRBL  (12) 

DIMENSION  VVRBL  (12),  FACTR  (12) 

DATA  (VVRBL(I),  1-1,12)/  / 

DATA  (VftBL(I),  I  -  1,12)/  / 

CALL  DRIVNF 
END 


This  program  sets  up  the  number  of  parameters,  the  parameter 
values,  the  independent  parameter  and  the  number  of  degrees-of-freedom  of 
the  problem.  The  variables  that  need  to  be  defined  are:- 

VVRBL  -  Array  containing  the  values  of  the  parameters  (allowed  up 
to  12  parameters). 

FACTR  -  Array  of  values  used  to  factorize  the  parameters  (set 
FACTR (I)  -  1/VVRBL(I ) ) . 


IVRBL  -  Index  of  the  Independent  parameter. 

IV1 ,IV2  -  Indloes  of  the  dependent  parameters. 

NM  -  Number  of  degrees-of-freedoa  of  the  system.  (NM  <  10). 

-  Charaoter  array  containing  a  descriptor  of  the 
corresponding  parameter,  up  to  10  characters  are  allowed 
(i.e.  VNBL(I)  contains  descriptor  of  parameter  VVNBL(I)). 


VRBL 
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Subroutine  SYSCAL 


COMMON/ VR/VVRBL,  FACTR,  IVRBL,  IV1 ,  IV2,  NM 

COMMON/SYS/SYSM 

DIMENSION  VVRBL(12),  FACTR (12) 

COMPLEX  SYSM(IO.IO) 


Given  the  array  of  parameters,  VVRBL,  this  routine  will  calculate  the 
system  matrix  S. 


SYSM 


Complex  array  which  contains  the  system  matrix  S 
order  of  matrix  S  is  NM  x  NM,  in  this  case  NM  <  10. 


The 


APPENDIX  C 

Two  Degree-of-Freedom  Example 


Table  2  gives  a  list  of  values  of  the  system  parameters  used  to  solve 
eqn  13.  Listings  of  the  user  supplied  routines  are  given  below. 


PROGRAM  PREPFLSO 

C 

C  This  sets  up  o  two  degree  of  freedom  flutter 

C  problem. (i. e.  eqn.14)  (see  DOUELL  et  ol.»  pBO) 

C 

C  This  program  starts  VAFLUTSQLV  by- 

C  (1)  specifying  the  no.  of  D.O.F.  of  the  system. 

C  (2)  specifying  the  parameters  of  the  system  and 

C  their  initial  values. 

C  (3)  specifying  the  initial  independent  parameter 

C 

COMMON/ VR/VVRBL,  FACTR*  IVRBL,  IV1,  IV2,  NM 

COMMON/LAB/VRBL 

CHARACTER*10  VRBL<12> 

DIMENSION  VVRBL(12>,FACTR<12> 

DATA  <VVRBL(I>, 1=1, 12)/ 

1  10.  ,0.5,  0  5,  0.20,  0.4,  1.0,0.0,0.5,6.28318,0.0, 

2  0.0,  0.  0/ 

DATA  <VRBL(I>, 1=1, 10)/ 

'  BETA  ','  FREQ. RAT.', 

'  RAD . GYRA  ' , '  X  MASS 
'  R  LIFT  ','  VELOCITY  ', 

'  DEC. RATIO','  FREQ 
'  DC1/DA  ','  FREQ  '/ 

DATA  PI/3.14159265/ 


c 

INDEX 

VARIABLE 

VALUE 

c 

1 

BETA 

10. 

c 

2 

FREQ. RATIO 

0.  5 

c 

3 

RAD.GYRAT. 

0.5 

c 

4 

MASS  MOM. ARM 

0.20 

c 

5 

LIFT  MOM. ARM 

0.  4 

c 

6 

NON. DIM. VEL. 

0.001 

c 

7 

NON. DIM. DECAY  RA. 

0.001 

c 

8 

NON. DIM. FREQ1 

0.5 

c 

9 

LIFT  CURVE  SL. 

2PI 

C**  D.O.F.  OF  THE  SYSTEM  ** 

NM=2 

C**  SET  UP  FACTORIZING  VALUES  ** 

DO  10  1=1, 12 
10  FACTR < I )*0. 0 

DO  20  1-1, 10 

IF( VVRBL< I ) . EQ. 0. 0)  THEN 

FACTR< I )»1 . 0 

ELSE 

FACTR ( I ) = 1 . /VVRBL ( I ) 

ENDIF 

20  CONTINUE 

Ct«  8ET  UP  INITIAL  INDEP  VARBLE  ** 
IVRBL *6 

CALL  DRIVNF 


STOP 

END 


SUBROUTINE  SYSCAL 

C 

C  This  routine  sets  up  the  system  matrix  S. 

C 

COMMON/VR/WRBL,  FACTR,  IVRBL,  I VI,  IV2,  NM 

COMMON/SYS/  SYSM 

DIMENSION  WRBL <  1 2 )  >  FACTR <  1 2 > 

COMPLEX  SYSM (10, 10) , M<2, 2) . S<2, 2) , A<2, 2> . CHI, 1 1 
DATA  PI/3. 14159265/ 

I I=CMPLX<0. 0, 1. 0> 

C**  MASS/INERTIA  MATRIX  ** 

M<1, 1>=1. 0 
M(l>  2>=0. 20 
M<2» 1>=M<1, 2) 

M(2> 2 ) =0 . 25 

CHI=WRBL<7)  +  <II*WRBL<8>> 

C**  STIFFNESS  MATRIX  ** 

S<1, 1>=1. 0 
S(l>  2)=0. 0 
S<2, 1  )=0. 0 

S<2,2)  =  <WRBL<3>**2.  >*U . /<WRBL<2)**2.  )) 

C**  AERO.  MATRIX  ** 

All, 1>=CHI 
A  ( 1  <  2 )  =VVRBL  <  6 ) 

A<2» 1>=-CHI*0. 4 
A<2,  2>=-0.  4*WRBL<6) 

C**  CALC.  SYSTEM  MATRIX  ** 

DO  10  1=1,2 
DO  10  J=1 , 2 

Mil, J>=31. 4#<CHI**2)*M<I, J> 

S<I, J)=7. 85*S< I , J) 

All,  J>»6.  28*A< I«  J)*<WRBL<6) ) 

10  SYSM( I, J)=M(I*  J)+S< I, J)+A< I,  J) 


RETURN 

END 


TABLE  1 


Summary  of  Flutter  Solver  Routines 


Subroutines  in 
VAFLUTSOLV  File 

j  Description 

INIT 

Input  data,  start  program 

OUTP 

Output  data 

CMAND 

Graphics  input/output  data 

DISP 

Graphics,  Output  data 

OLABEL 

Graphics,  Output  data 

TPTEXT 

Graphics,  Output  data 

DRIVNF 

Main  solving  routine 

VORT 

Called  by  DRIVNF 

ANG 

Called  by  DRIVNF 

CIRCLE 

Called  by  DRIVNF 

NF 

Called  by  DRIVNF,  calls  F2V 

FG 

Calls  SYSCAL  to  determine  S 

F2V 

Refines  estimates  x^,  xk 

CMINV 

Input  A  determine  A  ,  |a| 

Other  Routines  | 

i 

1 

SYSCAL  | 

1 

|  User  supplied  routine  to  determine  S 

PREPFLSQ  | 

|  User  supplied  routine  to  set  up  problem 

(MOTE:  These  files  are  contained  in  sub-directory  FLUTSOLVER) 


List  of  Parameter 


1  Parameter 

Description 

(VRBL) 

0 

| 

Mass  Parameter 

Uh  /(I) 

n  a 

Freq.  Ratio 

r 

i  a 

Radius  of  Gyration 

X 

a 

1 

Mass  Moment  Arm 

\ 

i  e 

Lift  Moment  Arm 

1  .  5 

Non-dim  Velocity 

*  {  w/aj 

t  a 

Non-dim.  Freq. 

,  Acj/w 

a 

Decay  Rate 

!  3CL/^a 

Lift  Curve  Slope 

2 


uea  used  to  Solve  EQN  (1*0 


Parameter 
Value  or 
Initial  Guess 
(VVRBL) 


Index 

of 

VVRBL 


Factorising 
Value 
(FACTR ) 


•  Solution 
x  Estimate 
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METHOD  FOR  PREDICTING  THE  NEW  INITIAL  ESTIMATE 
OF  (SIMILAR  PROCEDURE  FOR  xR) 
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BLOCK  DIAGRAM  OF  THE  SUBROUTINES 
USED  TO  SOLVE  THE  FLUTTER  EQUATIONS 
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PIG.  4<a)  VARIATION  OF  NON-DIMENSIONAL  FREQUENCY  AND  DAMPING 
WITH  NON-DIMENSIONAL  VELOCITY  FOR  THE  FIRST 
NATURAL  MODE 
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An  Interactive  direct  flutter  solving  routine  has  been 
installed  on  ARL's  ELXSI  computer.  Interactive  graphics  routines 
using  DI-3000  graphics,  software,  have  been  Incorporated  to  give 
the  user  a  progressive  picture  of  the  solution.  Estimates  of 
subcrltlcal  response  data  may  also  be  determined  using  this 
method.  A  description  of  the  method  and  the  associated  software 
is  presented  here.  Also  included  is  a  two-dimensional  flutter 
problem,  using  quasi-steady  aerodynamics,  which  is  solved  by  this 
direct  flutter  solver.  Yty <*>,■>, Js  1 


